# Copyright 2023 The GPU4PySCF Authors. All Rights Reserved.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

'''
SMD solvent model, copied from GPU4PYSCF with modification for CPU
'''

import numpy as np
from pyscf import lib, gto
from pyscf.data import radii
from pyscf.dft import gen_grid
from pyscf.solvent import pcm
from pyscf.solvent._attach_solvent import _for_scf
from pyscf.lib import logger

@lib.with_doc(_for_scf.__doc__)
def smd_for_scf(mf, solvent_obj=None, dm=None):
    if solvent_obj is None:
        solvent_obj = SMD(mf.mol)
    return _for_scf(mf, solvent_obj, dm)

# Inject PCM to SCF, TODO: add it to other methods later
from pyscf import scf
scf.hf.RHF.SMD = smd_for_scf
scf.uhf.UHF.SMD = smd_for_scf
hartree2kcal = 627.509451

# database from https://comp.chem.umn.edu/solvation/mnsddb.pdf
# solvent name: [n, n25, alpha, beta, gamma, epsilon, phi, psi)
solvent_db = {
    '1,1,1-trichloroethane':[1.4379, 1.4313, 0.0, 0.09, 36.24, 7.0826, 0.0, 0.60],
    '1,1,2-trichloroethane':[1.4717, 1.4689, 0.13, 0.13, 48.97, 7.1937, 0.0, 0.60],
    '1,2,4-trimethylbenzene':[1.5048, 1.5024, 0.0, 0.19, 42.03, 2.3653, 0.667, 0.0],
    '1,2-dibromoethane':[1.5387, 1.5364, 0.10, 0.17, 56.93, 4.9313, 0.0, 0.5],
    '1,2-dichloroethane':[1.4448, 1.4425, 0.10, 0.11, 45.86, 10.125, 0.0, 0.5],
    '1,2-ethanediol':[1.4318, 1.4306, 0.58, 0.78, 69.07, 40.245, 0.0, 0.5],
    '1,4-dioxane':[1.4224, 1.4204, 0.00, 0.64, 47.14, 2.2099, 0.0, 0.0],
    '1-bromo-2-methylpropane':[1.4348, 1.4349, 0.00, 0.12, 34.69, 7.7792, 0.0, 0.2],
    '1-bromooctane':[1.4524, 1.4500, 0.0, 0.12, 41.28, 5.0244, 0.0, 0.111],
    '1-bromopentane':[1.4447, 1.4420, 0.00, 0.12, 38.7, 6.269, 0.0, 0.167],
    '1-bromopropane':[1.4343, 1.4315, 0.0, 0.12, 36.36, 8.0496, 0.0, 0.250],
    '1-butanol':[1.3993, 1.3971, 0.37, 0.48, 35.88, 17.332, 0.0, 0.0],
    '1-chlorohexane':[1.4199, -1, 0.0, 0.10, 37.03, 5.9491, 0.0, 0.143],
    '1-chloropentane':[1.4127, 1.4104, 0.0, 0.1, 35.12, 6.5022, 0.0, 0.167],
    '1-chloropropane':[1.3879, 1.3851, 0.0, 0.1, 30.66, 8.3548, 0.0, 0.25],
    '1-decanol':[1.4372, 1.4353, 0.37, 0.48, 41.04, 7.5305, 0.0, 0.0],
    '1-fluorooctane':[1.3935, 1.3927, 0.0, 0.10, 33.92, 3.89, 0.0, 0.111],
    '1-heptanol':[1.4249, 1.4224, 0.37, 0.48, 38.5, 11.321, 0.0, 0.0],
    '1-hexanol':[1.4178, 1.4162, 0.37, 0.48, 37.15, 12.51, 0.0, 0.0],
    '1-hexene':[1.3837, 1.385, 0.00, 0.07, 25.76, 2.0717, 0.0, 0.0],
    '1-hexyne':[1.3989, 1.3957, 0.12, 0.10, 28.79, 2.615, 0.0, 0.0],
    '1-iodobutane':[1.5001, 1.4958, 0.00, 0.15, 40.65, 6.173, 0.0, 0.0],
    '1-iodohexadecane':[1.4806, -1, 0.00, 0.15, 46.48, 3.5338, 0.0, 0.0],
    '1-iodopentane':[1.4959, -1, 0.00, 0.15, 41.56, 5.6973, 0.0, 0.0],
    '1-iodopropane':[1.5058, 1.5027, 0.00, 0.15, 41.45, 6.9626, 0.0, 0.0],
    '1-nitropropane':[1.4018, 1.3996, 0.00, 0.31, 43.32, 23.73, 0.0, 0.0],
    '1-nonanol':[1.4333, 1.4319, 0.37, 0.48, 40.14, 8.5991, 0.0, 0.0],
    '1-octanol':[1.4295, 1.4279, 0.37, 0.48, 39.01, 9.8629, 0.0, 0.0],
    '1-pentanol':[1.4101, 1.4080, 0.37, 0.48, 36.5, 15.13, 0.0, 0.0],
    '1-pentene':[1.3715, 1.3684, 0.00, 0.07, 22.24, 1.9905, 0.0, 0.0],
    '1-propanol':[1.3850, 1.3837, 0.37, 0.48, 33.57, 20.524, 0.0, 0.0],
    '2,2,2-trifluoroethanol':[1.2907, -1, 0.57, 0.25, 42.02, 26.726, 0.0, 0.5],
    '2,2,4-trimethylpentane':[1.3915, 1.3889, 0.00, 0.00, 26.38, 1.9358, 0.0, 0.0],
    '2,4-dimethylpentane':[1.3815, 1.3788, 0.0, 0.00, 25.42, 1.8939, 0.0, 0.0],
    '2,4-dimethylpyridine':[1.5010, 1.4985, 0.0, 0.63, 46.86, 9.4176, 0.625, 0.0],
    '2,6-dimethylpyridine':[1.4953, 1.4952, 0.0, 0.63, 44.64, 7.1735, 0.625, 0.0],
    '2-bromopropane':[1.4251, 1.4219, 0.00, 0.14, 33.46, 9.3610, 0.0, 0.25],
    '2-butanol':[1.3978, 1.3949, 0.33, 0.56, 32.44, 15.944, 0.0, 0.0],
    '2-chlorobutane':[1.3971, 1.3941, 0.00, 0.12, 31.1, 8.3930, 0.0, 0.2],
    '2-heptanone':[1.4088, 1.4073, 0.0, 0.51, 37.6, 11.658, 0.0, 0.0],
    '2-hexanone':[1.4007, 1.3987, 0.0, 0.51, 36.63, 14.136, 0.0, 0.0],
    '2-methoxyethanol':[1.4024, 1.4003, 0.30, 0.84, 44.39, 17.2, 0.0, 0.0],
    '2-methyl-1-propanol':[1.3955, 1.3938, 0.37, 0.48, 32.38, 16.777, 0.0, 0.0],
    '2-methyl-2-propanol':[1.3878, 1.3852, 0.31, 0.60, 28.73, 12.47, 0.0, 0.0],
    '2-methylpentane':[1.3715, 1.3687, 0.0, 0.00, 24.3, 1.89, 0.0, 0.0],
    '2-methylpyridine':[1.4957, 1.4984, 0.0, 0.58, 47.5, 9.9533, 0.714, 0.0],
    '2-nitropropane':[1.3944, 1.3923, 0.00, 0.33, 42.16, 25.654, 0.00, 0.0],
    '2-octanone':[1.4151, 1.4133, 0.00, 0.51, 37.29, 9.4678, 0.0, 0.0],
    '2-pentanone':[1.3895, 1.3885, 0.00, 0.51, 33.46, 15.200, 0.0, 0.0],
    '2-propanol':[1.3776, 1.3752, 0.33, 0.56, 30.13, 19.264, 0.0, 0.0],
    '2-propen-1-ol':[1.4135, 0.00, 0.38, 0.48, 36.39, 19.011, 0.0, 0.0],
    'E-2-pentene':[1.3793, 1.3761, 0.00, 0.07, 23.62, 2.051, 0.0, 0.0],
    '3-methylpyridine':[1.5040, 1.5043, 0.0, 0.54, 49.61, 11.645, 0.714, 0.0],
    '3-pentanone':[1.3924, 1.3905, 0.00, 0.51, 35.61, 16.78, 0.0, 0.0],
    '4-heptanone':[1.4069, 1.4045, 0.00, 0.51, 35.98, 12.257, 0.0, 0.0],
    '4-methyl-2-pentanone':[1.3962, 1.394, 0.0, 0.51, 33.83, 12.887, 0.0, 0.0],
    '4-methylpyridine':[1.5037, 1.503, 0.0, 0.54, 50.17, 11.957, 0.714, 0.0],
    '5-nonanone':[1.4195, 0.00, 0.0, 0.51, 37.83, 10.6, 0.0, 0.0],
    'acetic acid':[1.3720, 1.3698, 0.61, 0.44, 39.01, 6.2528, 0.0, 0.0],
    'acetone':[1.3588, 1.3559, 0.04, 0.49, 33.77, 20.493, 0.0, 0.0],
    'acetonitrile':[1.3442, 1.3416, 0.07, 0.32, 41.25, 35.688, 0.0, 0.0],
    'acetophenone':[1.5372, 1.5321, 0.00, 0.48, 56.19, 17.44, 0.667, 0.0],
    'aniline':[1.5863, 1.5834, 0.26, 0.41, 60.62, 6.8882, 0.857, 0.0],
    'anisole':[1.5174, 1.5143, 0.0, 0.29, 50.52, 4.2247, 0.75, 0.0],
    'benzaldehyde':[1.5463, 1.5433, 0.0, 0.39, 54.69, 18.220, 0.857, 0.0],
    'benzene':[1.5011, 1.4972, 0.0, 0.14, 40.62, 2.2706, 1.0, 0.0],
    'benzonitrile':[1.5289, 1.5257, 0.0, 0.33, 55.83, 25.592, 0.75, 0.0],
    'benzylalcohol':[1.5396, 1.5384, 0.33, 0.56, 52.96, 12.457, 0.75, 0.0],
    'bromobenzene':[1.5597, 1.5576, 0.0, 0.09, 50.72, 5.3954, 0.857, 0.143],
    'bromoethane':[1.4239, 1.4187, 0.0, 0.12, 34.0, 9.01, 0.0, 0.333],
    'bromoform':[1.6005, 1.5956, 0.15, 0.06, 64.58, 4.2488, 0.0, 0.75],
    'butanal':[1.3843, 1.3766, 0.0, 0.45, 35.06, 13.45, 0.0, 0.0],
    'butanoic acid':[1.3980, 1.3958, 0.60, 0.45, 37.49, 2.9931, 0.0, 0.0],
    'butanone':[1.3788, 1.3764, 0.00, 0.51, 34.5, 18.246, 0.0, 0.0],
    'butanonitrile':[1.3842, 1.382, 0.0, 0.36, 38.75, 24.291, 0.0, 0.0],
    'butylethanoate':[1.3941, 1.3923, 0.0, 0.45, 35.81, 4.9941, 0.0, 0.0],
    'butylamine':[1.4031, 1.3987, 0.16, 0.61, 33.74, 4.6178, 0.0, 0.0],
    'n-butylbenzene':[1.4898, 1.4874, 0.0, 0.15, 41.33, 2.36, 0.6, 0.0],
    'sec-butylbenzene':[1.4895, 1.4878, 0.0, 0.16, 40.35, 2.3446, 0.60, 0.0],
    'tert-butylbenzene':[1.4927, 1.4902, 0.0, 0.16, 39.78, 2.3447, 0.6, 0.0],
    'carbon disulfide':[1.6319, 1.6241, 0.0, 0.07, 45.45, 2.6105, 0.0, 0.0],
    'carbon tetrachloride':[1.4601, 1.4574, 0.00, 0.00, 38.04, 2.2280, 0.0, 0.8],
    'chlorobenzene':[1.5241, 1.5221, 0.0, 0.07, 47.48, 5.6968, 0.857, 0.143],
    'chloroform':[1.4459, 1.4431, 0.15, 0.02, 38.39, 4.7113, 0.0, 0.75],
    'a-chlorotoluene':[1.5391, 0.0, 0.0, 0.33, 53.04, 6.7175, 0.75, 0.125],
    'o-chlorotoluene':[1.5268, 1.5233, 0.00, 0.07, 47.43, 4.6331, 0.75, 0.125],
    'm-cresol':[1.5438, 1.5394, 0.57, 0.34, 51.37, 12.44, 0.75, 0.0],
    'o-cresol':[1.5361, 1.5399, 0.52, 0.30, 53.11, 6.76, 0.75, 0.0],
    'cyclohexane':[1.4266, 1.4235, 0.00, 0.00, 35.48, 2.0165, 0.0, 0.0],
    'cyclohexanone':[1.4507, 1.4507, 0.00, 0.56, 49.76, 15.619, 0.00, 0.0],
    'cyclopentane':[1.4065, 1.4036, 0.00, 0.00, 31.49, 1.9608, 0.0, 0.0],
    'cyclopentanol':[1.4530, -1, 0.32, 0.56, 46.8, 16.989, 0.0, 0.0],
    'cyclopentanone':[1.4366, 1.4347, 0.00, 0.52, 47.21, 13.58, 0.0, 0.0],
    'decalin (cis/trans mixture)':[1.4753, 1.472, 0.00, 0.00, 43.82, 2.196, 0.0, 0.0],
    'cis-decalin':[1.4810, 1.4788, 0.00, 0.00, 45.45, 2.2139, 0.0, 0.0],
    'n-decane':[1.4102, 1.4094, 0.00, 0.00, 33.64, 1.9846, 0.0, 0.0],
    'dibromomethane':[1.5420, 1.5389, 0.10, 0.10, 56.21, 7.2273, 0.0, 0.667],
    'butylether':[1.3992, 1.3968, 0.00, 0.45, 35.98, 3.0473, 0.0, 0.0],
    'o-dichlorobenzene':[1.5515, 1.5491, 0.00, 0.04, 52.72, 9.9949, 0.75, 0.25],
    'E-1,2-dichloroethene':[1.4454, 1.4435, 0.09, 0.05, 37.13, 2.14, 0.0, 0.5],
    'Z-1,2-dichloroethene':[1.4490, 1.4461, 0.11, 0.05, 39.8, 9.2, 0.0, 0.5],
    'dichloromethane':[1.4242, 1.4212, 0.10, 0.05, 39.15, 8.93, 0.0, 0.667],
    'diethylether':[1.3526, 1.3496, 0.00, 0.41, 23.96, 4.2400, 0.0, 0.0],
    'diethylsulfide':[1.4430, 1.4401, 0.00, 0.32, 35.36, 5.723, 0.0, 0.0],
    'diethylamine':[1.3864, 1.3825, 0.08, 0.69, 28.57, 3.5766, 0.0, 0.0],
    'diiodomethane':[1.7425, 1.738, 0.05, 0.23, 95.25, 5.32, 0.0, 0.0],
    'diisopropyl ether':[1.3679, 1.3653, 0.00, 0.41, 24.86, 3.38, 0.0, 0.0],
    'cis-1,2-dimethylcyclohexane':[1.4360, 1.4336, 0.00, 0.00, 36.28, 2.06, 0.0, 0.0],
    'dimethyldisulfide':[1.5289, 1.522, 0.00, 0.28, 48.06, 9.6, 0.0, 0.0],
    'N,N-dimethylacetamide':[1.4380, 1.4358, 0.00, 0.78, 47.62, 37.781, 0.0, 0.0],
    'N,N-dimethylformamide':[1.4305, 1.4280, 0.00, 0.74, 49.56, 37.219, 0.0, 0.0],
    'dimethylsulfoxide':[1.4783, 1.4783, 0.00, 0.88, 61.78, 46.826, 0.0, 0.0],
    'diphenylether':[1.5787, -1, 0.00, 0.20, 38.5, 3.73, 0.923, 0.0],
    'dipropylamine':[1.4050, 1.4018, 0.08, 0.69, 32.11, 2.9112, 0.0, 0.0],
    'n-dodecane':[1.4216, 1.4151, 0.00, 0.00, 35.85, 2.0060, 0.0, 0.0],
    'ethanethiol':[1.4310, 1.4278, 0.00, 0.24, 33.22, 6.667, 0.0, 0.0],
    'ethanol':[1.3611, 1.3593, 0.37, 0.48, 31.62, 24.852, 0.0, 0.0],
    'ethylethanoate':[1.3723, 1.3704, 0.00, 0.45, 33.67, 5.9867, 0.0, 0.0],
    'ethylmethanoate':[1.3599, 1.3575, 0.00, 0.38, 33.36, 8.3310, 0.0, 0.0],
    'ethylphenylether':[1.5076, 1.5254, 0.00, 0.32, 46.65, 4.1797, 0.667, 0.0],
    'ethylbenzene':[1.4959, 1.4932, 0.00, 0.15, 41.38, 2.4339, 0.75, 0.0],
    'fluorobenzene':[1.4684, 1.4629, 0.00, 0.10, 38.37, 5.42, 0.857, 0.143],
    'formamide':[1.4472, 1.4468, 0.62, 0.60, 82.08, 108.94, 0.0, 0.0],
    'formicacid':[1.3714, 1.3693, 0.75, 0.38, 53.44, 51.1, 0.0, 0.0],
    'n-heptane':[1.3878, 1.3855, 0.00, 0.00, 28.28, 1.9113, 0.0, 0.0],
    'n-hexadecane':[1.4345, 1.4325, 0.00, 0.00, 38.93, 2.0402, 0.0, 0.0],
    'n-hexane':[1.3749, 1.3722, 0.00, 0.00, 25.75, 1.8819, 0.0, 0.0],
    'hexanoicacid':[1.4163, 1.4146, 0.60, 0.45, 39.65, 2.6, 0.0, 0.0],
    'iodobenzene':[1.6200, 1.6172, 0.00, 0.12, 55.72, 4.5470, 0.857, 0.0],
    'iodoethane':[1.5133, 1.5100, 0.00, 0.15, 40.96, 7.6177, 0.0, 0.0],
    'iodomethane':[1.5380, 1.5270, 0.00, 0.13, 43.67, 6.8650, 0.0, 0.0],
    'isopropylbenzene':[1.4915, 1.4889, 0.00, 0.16, 39.85, 2.3712, 0.667, 0.0],
    'p-isopropyltoluene':[1.4909, 1.4885, 0.00, 0.19, 38.34, 2.2322, 0.600, 0.0],
    'mesitylene':[1.4994, 1.4968, 0.00, 0.19, 39.65, 2.2650, 0.667, 0.0],
    'methanol':[1.3288, 1.3265, 0.43, 0.47, 31.77, 32.613, 0.0, 0.0],
    'methylbenzoate':[1.5164, 1.5146, 0.00, 0.46, 53.5, 6.7367, 0.600, 0.0],
    'methylbutanoate':[1.3878, 1.3847, 0.00, 0.45, 35.44, 5.5607, 0.0, 0.0],
    'methylethanoate':[1.3614, 1.3589, 0.00, 0.45, 35.59, 6.8615, 0.0, 0.0],
    'methylmethanoate':[1.3433, 1.3415, 0.00, 0.38, 35.06, 8.8377, 0.0, 0.0],
    'methylpropanoate':[1.3775, 1.3742, 0.00, 0.45, 35.18, 6.0777, 0.0, 0.0],
    'N-methylaniline':[1.5684, 1.5681, 0.17, 0.43, 53.11, 5.9600, 0.75, 0.0],
    'methylcyclohexane':[1.4231, 1.4206, 0.00, 0.00, 33.52, 2.024, 0.0, 0.0],
    'N-methylformamide(E/Zmixture)':[1.4319, 1.4310, 0.40, 0.55, 55.44, 181.56, 0.0, 0.0],
    'nitrobenzene':[1.5562, 1.5030, 0.00, 0.28, 57.54, 34.809, 0.667, 0.0],
    'nitroethane':[1.3917, 1.3897, 0.02, 0.33, 46.25, 28.29, 0.0, 0.0],
    'nitromethane':[1.3817, 1.3796, 0.06, 0.31, 52.58, 36.562, 0.0, 0.0],
    'o-nitrotoluene':[1.5450, 1.5474, 0.0, 0.27, 59.12, 25.669, 0.6, 0.0],
    'n-nonane':[1.4054, 1.4031, 0.0, 0.0, 32.21, 1.9605, 0.0, 0.0],
    'n-octane':[1.3974, 1.3953, 0.0, 0.0, 30.43, 1.9406, 0.0, 0.0],
    'n-pentadecane':[1.4315, 1.4298, 0.0, 0.0, 38.34, 2.0333, 0.0, 0.0],
    'pentanal':[1.3944, 1.3917, 0.0, 0.4, 36.62, 10.0, 0.0, 0.0],
    'n-pentane':[1.3575, 1.3547, 0.0, 0.0, 22.3, 1.8371, 0.0, 0.0],
    'pentanoic acid':[1.4085, 1.4060, 0.60, 0.45, 38.4, 2.6924, 0.0, 0.0],
    'pentyl ethanoate':[1.4023, -1.0, 0.0, 0.45, 36.23, 4.7297, 0.0, 0.0],
    'pentylamine':[1.448, 1.4093, 0.16, 0.61, 35.54, 4.2010, 0.0, 0.0],
    'perfluorobenzene':[1.3777, 1.3761, 0.00, 0.00, 31.74, 2.029, 0.5, 0.5],
    'propanal':[1.3636, 1.3593, 0.00, 0.45, 32.48, 18.5, 0.0, 0.0],
    'propanoic acid':[1.3869, 1.3848, 0.60, 0.45, 37.71, 3.44, 0.0, 0.0],
    'propanonitrile':[1.3655, 1.3633, 0.02, 0.36, 38.5, 29.324, 0.0, 0.0],
    'propyl ethanoate':[1.3842, 1.3822, 0.0, 0.45, 34.26, 5.5205, 0.0, 0.0],
    'propylamine':[1.3870, 1.3851, 0.16, 0.61, 31.31, 4.9912, 0.0, 0.0],
    'pyridine':[1.5095, 1.5073, 0.0, 0.52, 52.62, 12.978, 0.833, 0.0],
    'tetrachloroethene':[1.5053, 1.5055, 0.0, 0.0, 45.19, 2.268, 0.0, 0.667],
    'tetrahydrofuran':[1.4050, 1.4044, 0.0, 0.48, 39.44, 7.4257, 0.0, 0.0],
    'tetrahydrothiophene-S,S-dioxide':[1.4833, -1.0, 0.0, 0.88, 87.49, 43.962, 0.0, 0.0],
    'tetralin':[1.5413, 1.5392, 0.0, 0.19, 47.74, 2.771, 0.6, 0.0],
    'thiophene':[1.5289, 1.5268, 0.0, 0.15, 44.16, 2.7270, 0.8, 0.0],
    'thiophenol':[1.5893, 1.580, 0.09, 0.16, 55.24, 4.2728, 0.857, 0.0],
    'toluene':[1.4961, 1.4936, 0.0, 0.14, 40.2, 2.3741, 0.857, 0.0],
    'trans-decalin':[1.4695, 1.4671, 0.0, 0.0, 42.19, 2.1781, 0.0, 0.0],
    'tributylphosphate':[1.4224, 1.4215, 0.0, 1.21, 27.55, 8.1781, 0.0, 0.0],
    'trichloroethene':[1.4773, 1.4556, 0.08, 0.03, 41.45, 3.422, 0.0, 0.6],
    'triethylamine':[1.4010, 1.3980, 0.0, 0.79, 29.1, 2.3832, 0.0, 0.0],
    'n-undecane':[1.4398, 1.4151, 0.0, 0.0, 34.85, 1.991, 0.0, 0.0],
    'water':[1.3328, 1.3323, 0.82, 0.35, -1.0, 78.355, -1.0, -1.0],
    'xylene (mixture)':[1.4995, 1.4969, 0.0, 0.16, 41.38, 2.3879, 0.75, 0.0],
    'm-xylene':[1.4972, 1.4946, 0.0, 0.16, 40.98, 2.3478, 0.75, 0.0],
    'o-xylene':[1.5055, 1.5029, 0.0, 0.16, 42.83, 2.5454, 0.75, 0.0],
    'p-xylene':[1.4958, 1.4932, 0.0, 0.16, 40.32, 2.2705, 0.75, 0.0],
    '':[0]*8
}

def smd_radii(alpha):
    '''
    eq. (16)
    use smd radii if defined
    use Bondi radii if defined
    use 2.0 otherwise
    '''
    radii_table = radii.VDW.copy() * radii.BOHR
    radii_table[1] = 1.20
    radii_table[6] = 1.85
    radii_table[7] = 1.89
    if alpha >= 0.43:
        r = 1.52
    else:
        r = 1.52 + 1.8 * (0.43 - alpha)
    radii_table[8] = r
    radii_table[9] = 1.73
    radii_table[14] = 2.47
    radii_table[15] = 2.12
    radii_table[16] = 2.49
    radii_table[17] = 2.38
    #radii_table[35] = 3.06 # original SMD
    # following value from SMD18
    # https://chemistry-europe.onlinelibrary.wiley.com/doi/10.1002/chem.201803652
    radii_table[35] = 2.60
    radii_table[53] = 2.74
    return radii_table/radii.BOHR

import ctypes
from pyscf.lib import load_library
try:
    libsolvent = load_library('libsolvent')
except (IOError, NameError):
    libsolvent = None

def get_cds_legacy(smdobj):
    if libsolvent is None:
        raise RuntimeError(
            'SMD module is not available. '
            'You can compile this module with cmake option "-DENABLE_SMD=ON"')

    mol = smdobj.mol
    natm = mol.natm
    soln, _, sola, solb, solg, _, solc, solh = smdobj.solvent_descriptors
    #symbols = [mol.atom_s(ia) for ia in range(mol.natm)]
    charges = np.asarray(mol.atom_charges(), dtype=np.int32, order='F')
    coords = np.asarray(mol.atom_coords(unit='B'), dtype=np.float64, order='C')
    icds = 1 if smdobj.solvent.upper() == 'WATER' else 2
    dcds = np.empty([natm,3])
    mnsol_interface =  libsolvent.mnsol_interface_

    double_ndptr = np.ctypeslib.ndpointer(dtype=np.float64)
    int_ndptr = np.ctypeslib.ndpointer(dtype=np.int32)
    double_ptr = ctypes.POINTER(ctypes.c_double)
    int_ptr = ctypes.POINTER(ctypes.c_int)

    mnsol_interface.argtypes = [
        double_ndptr, int_ndptr,
        int_ptr,
        double_ptr, double_ptr, double_ptr, double_ptr, double_ptr, double_ptr,
        int_ptr,
        double_ptr, double_ptr, double_ndptr]
    natm = ctypes.byref(ctypes.c_int(natm))
    icds = ctypes.byref(ctypes.c_int(icds))
    soln = ctypes.byref(ctypes.c_double(soln))
    sola = ctypes.byref(ctypes.c_double(sola))
    solb = ctypes.byref(ctypes.c_double(solb))
    solg = ctypes.byref(ctypes.c_double(solg))
    solc = ctypes.byref(ctypes.c_double(solc))
    solh = ctypes.byref(ctypes.c_double(solh))
    gcds = ctypes.c_double()
    areacds = ctypes.c_double()

    mnsol_interface(coords, charges,
                    natm,
                    sola, solb, solc, solg, solh, soln,
                    icds,
                    ctypes.byref(gcds), ctypes.byref(areacds), dcds)
    return gcds.value / hartree2kcal, dcds

# Note: in various places, SMD instance is not explictly tested. It is checked
# by the statement "isinstance(solvent, PCM)"
class SMD(pcm.PCM):
    '''
    SMD Solvent Model

    Attributes:
    ----------
    method : str
        No effects. It is set to 'SMD' as a placeholder

    vdw_scale : float
        A scaling factor for van der Waals radii. Default is 1.0.

    r_probe : float
        An additional radius (in Angstrom) added to the van der Waals radii.
        Default is 0.4 Angstrom.

    radii_table : dict
        Custom van der Waals radii for each element. By default, scaled van der Waals radii
        from `vdw_scale` and `r_probe` are used.

    sasa_ng : int
        The number of quadrature grids used for calculating the Solvent
        Accessible Surface Area (SASA). Default is 590.

    solvent : str
        The name of the solvent, which is used to determine the dielectric constant and other
        relevant parameters. Supported solvents can be accessed via the variable
        `pyscf.solvent.smd.solvent_db`.

    frozen : bool
        Whether to freeze the potential produced by the solvent during SCF iterations or
        other convergence processes. When frozen=True is set, the solvent is
        assumed to respond slowly, while the electron density relaxes quickly.
        Default is False.

    max_cycle : int
        The maximum number of iterations to relax the solvent.

    conv_tol : float
        The convergence tolerance for total energy during solvent relaxation.

    equilibrium_solvation : bool
        Affects TDDFT and other excited state computations. Controls whether the solvent
        relaxes rapidly with respect to the electron density of the excited state.
        For vertical excitations, it is recommended to set this to False, as the solvent
        typically does not fully relax. In some software packages (e.g., Q-Chem),
        non-equilibrium solvation is applied with an optical dielectric constant of
        eps=1.78. Default is False.

    state_id : int
        Specifies the target state in excited state calculations.
        `state_id=0` corresponds to the ground state, while `state_id=1` corresponds
        to the first excited state. Default is 0.

    Saved Results:
    --------------
    e_cds : float
        Cavitation, Dispersion and Solvent energy

    Intermediate Attributes:
    ------------------------
    These attributes are generated during calculations and should not be modified.
    Additionally, they may not be compatible between GPU and CPU implementations.

    - surface
    - _intermediates
    - v_grids_n
    - solvent_descriptors
    '''

    _keys = {
        'method', 'e_cds', 'solvent_descriptors', 'r_probe', 'sasa_ng'
    }
    def __init__(self, mol, solvent=''):
        super().__init__(mol)
        self.vdw_scale = 1.0
        self.sasa_ng = 590 # quadrature grids for calculating SASA
        self.r_probe = 0.4/radii.BOHR
        self.method = 'SMD'  # use IEFPCM for electrostatic
        if solvent not in solvent_db:
            raise RuntimeError(f'{solvent} is not available in SMD')
        self._solvent = solvent
        self.solvent_descriptors = solvent_db[solvent]
        self.radii_table = None
        self.e_cds = None

    def build(self, ng=None):
        if self.radii_table is None:
            radii_table = smd_radii(self.solvent_descriptors[2])
        else:
            radii_table = self.radii_table
        logger.info(self, 'radii_table %s', radii_table*radii.BOHR)
        mol = self.mol
        if ng is None:
            ng = gen_grid.LEBEDEV_ORDER[self.lebedev_order]

        self.surface = pcm.gen_surface(mol, rad=radii_table, ng=ng)
        self._intermediates = {}
        F, A = pcm.get_F_A(self.surface)
        D, S = pcm.get_D_S(self.surface, with_S=True, with_D=True)

        epsilon = self.eps
        f_epsilon = (epsilon - 1.0)/(epsilon + 1.0)
        DA = D*A
        DAS = np.dot(DA, S)
        K = S - f_epsilon/(2.0*np.pi) * DAS
        R = -f_epsilon * (np.eye(K.shape[0]) - 1.0/(2.0*np.pi)*DA)
        intermediates = {
            'S': S,
            'D': D,
            'A': A,
            'K': K,
            'R': R,
            'f_epsilon': f_epsilon
        }
        self._intermediates.update(intermediates)

        charge_exp  = self.surface['charge_exp']
        grid_coords = self.surface['grid_coords']
        atom_coords = mol.atom_coords(unit='B')
        atom_charges = mol.atom_charges()

        int2c2e = mol._add_suffix('int2c2e')
        fakemol = gto.fakemol_for_charges(grid_coords, expnt=charge_exp**2)
        fakemol_nuc = gto.fakemol_for_charges(atom_coords)
        v_ng = gto.mole.intor_cross(int2c2e, fakemol_nuc, fakemol)
        self.v_grids_n = np.dot(atom_charges, v_ng)
        return self

    @property
    def solvent(self):
        return self._solvent

    @solvent.setter
    def solvent(self, solvent):
        self._solvent = solvent
        self.solvent_descriptors = solvent_db[solvent]
        self.radii_table = smd_radii(self.solvent_descriptors[2])
        self.eps = self.solvent_descriptors[5]
        self.reset()

    @property
    def sol_desc(self):
        return self.solvent_descriptors

    @sol_desc.setter
    def sol_desc(self, values):
        '''
        format of sol desc
        [n, n25, alpha, beta, gamma, epsilon, phi, psi]
        '''
        assert len(values) == 8
        self.solvent_descriptors = values
        self.radii_table = smd_radii(self.solvent_descriptors[2])
        self.eps = values[5]
        self.reset()

    def dump_flags(self, verbose=None):
        n, _, alpha, beta, gamma, _, phi, psi = self.solvent_descriptors
        logger.info(self, '******** %s ********', self.__class__)
        logger.info(self, 'lebedev_order = %s (%d grids per sphere)',
                    self.lebedev_order, gen_grid.LEBEDEV_ORDER[self.lebedev_order])
        logger.info(self, 'eps = %s'          , self.eps)
        logger.info(self, 'frozen = %s'       , self.frozen)
        logger.info(self, '---------- SMD solvent descriptors -------')
        logger.info(self, f'n     = {n}')
        logger.info(self, f'alpha = {alpha}')
        logger.info(self, f'beta  = {beta}')
        logger.info(self, f'gamma = {gamma}')
        logger.info(self, f'phi   = {phi}')
        logger.info(self, f'psi   = {psi}')
        logger.info(self, '--------------------- end ----------------')
        logger.info(self, 'equilibrium_solvation = %s', self.equilibrium_solvation)
        return self

    def get_cds(self):
        if self.e_cds is None:
            self.e_cds = get_cds_legacy(self)[0]
        return self.e_cds

    def nuc_grad_method(self, grad_method):
        raise DeprecationWarning('Use the make_grad_object function from '
                                 'pyscf.solvent.grad.smd instead.')

    def grad(self, dm, verbose=None):
        '''Computes the Jacobian for the energy associated with the solvent,
        including the derivatives of the solvent itsself and the interactions
        between the solvent and the charge density of the solute.
        '''
        from pyscf.solvent.grad.pcm import grad_qv, grad_nuc
        from pyscf.solvent.grad.smd import grad_solver, get_cds
        de_solvent = grad_qv(self, dm)
        de_solvent+= grad_solver(self, dm)
        de_solvent+= grad_nuc(self, dm)
        #de_cds     = get_cds(self.base.with_solvent)
        de_cds     = get_cds_legacy(self)[1]
        logger.info(self, 'Cavitation, Dispersion and Solvent structure contribution %s', de_cds)
        return de_solvent + de_cds

    def Hessian(self, hess_method):
        raise DeprecationWarning('Use the make_hess_object function from '
                                 'pyscf.solvent.hessian.smd instead.')

    def hess(self, dm):
        from pyscf.solvent.hessian.pcm import (
            analytical_hess_nuc, analytical_hess_qv, analytical_hess_solver)
        from pyscf.solvent.hessian.smd import get_cds
        de_solvent  =    analytical_hess_nuc(self, dm, verbose=self.verbose)
        de_solvent +=     analytical_hess_qv(self, dm, verbose=self.verbose)
        de_solvent += analytical_hess_solver(self, dm, verbose=self.verbose)
        de_cds = get_cds(self)
        logger.info(self, 'Cavitation, Dispersion and Solvent structure contribution %s', de_cds)
        return de_solvent + de_cds

    def reset(self, mol=None):
        super().reset(mol)
        self.e_cds = None
        return self
